Universal Model of Finite-Reynolds Number Turbulent Flow in Channels and Pipes 
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In this Letter we suggest a simple and physically transparent analytical model of pressure driven 
turbulent wall-bounded flows at high but finite Reynolds numbers Re. The model provides an 
accurate quantitative description of the profiles of the mean-velocity and Reynolds-stresses (second 
order correlations of velocity fluctuations) throughout the entire channel or pipe, for a wide range 
of Re, using only three Re- independent parameters. The model sheds light on the long-standing 
controversy between supporters of the century-old log-law theory of von-Karman and Prandtl and 
proposers of a newer theory promoting power laws to describe the intermediate region of the mean 
velocity profile. 



An important challenge in wall-bounded Newtonian 
turbulence is the description of the profiles of the 
mean velocity and second order correlation functions 
of turbulent-velocity fluctuations throughout the entire 
channel or pipe at relatively high but finite Reynolds 
numbers. To understand the issue, focus on a channel 
of width 2L between its parallel walls, where the incom- 
pressible fluid velocity C/(r, t) is decomposed into its av- 
erage (over time) and a fluctuating part 

U{r, t) = V{r) + u{r, t) , V{r) = {U{r, t)) 

. In a stationary plane channel flow with a constant pres- 
sure gradient p' = —dp/dx the only component of the 
mean velocity V is the stream-wise component Vx = V 
that depends on wall normal direction z only. Near the 
wall the mean velocity profiles for different Reynolds 
numbers exhibit data collapse once presented in wall 
units, where the Reynolds number Re,-, the normalized 
distance from the wall z'^ and the normalized mean ve- 
locity V'^{z~^) are deflned (for channels) by 



RCr 



= zRCr/L , V+ = VI y 



The classical theory of Prandtl and von-Karman for 
infinitely large Re,- is based on the assumption that the 
single characteristic scale in the problem is proportional 
to the distance from the (nearest) wall. It leads to the 
celebrated von-Karman log-law [1] 



y+(z+) = K-Mn(z+) + S, 



(1) 



which serves as a basis for the parametrization of turbu- 
lent flows near a wall in many engineering applications. 
On the face of it this law agrees with the data (see, e.g. 
Fig. 1) for relatively large 2+, say for z+ > 100, giving 
K ~ 0.4 and B ^ 5. The range of validly of the log-law 
is definitely restricted by the requirement C <C 1, where 
C = z/L (channel) ore ( = r/R (Pipe of radius R). For 
^ ~ 1 the global geometry becomes important leading 
to unavoidable deviations of V~^{Q from the log-law (1), 
known as the wake. 

The problem is that for finite Rct the corrections to 
the log-law (1) are in powers oi e = l/lnRe,- [5] and 



definitely cannot be neglected for the currently largest 
available direct numerical simulation (DNS) of channel 
flows (Re^ = 2003 [3] , giving e « 0.13). Even for Re^ 
approaching 500, 000 as in the Princeton Superpipe ex- 
periment [4], e « 0.08. This opens a Pandora box with 
various possibilities to revise the log-law (1) and to re- 
place it, as was suggested in [5], by a power law 



V+{z+) = C(Re^)(z+) 



(2) 



Here both C(ReT-) and 7(ReT-) were represented as 
asymptotic series expansions in e. The relative com- 
plexity of this proposition compared to the simplicity of 
Eq. (1) resulted in a less than enthusiastic response in the 
fluid mechanics community [6] , leading to a rather fierce 
controversy between the log-law camp and the power-law 
camp. Various attempts [4-9] to validate the log-law (1) 
or the alternative power-law (2) were based on extensive 
analysis of experimental data used to fit the velocity pro- 
files as a formal expansion in inverse powers of e or as 
composite expansions in both z^ and C. Note however 
that in the excellent fits presented, say in [9], one uses 
four adjustable parameters for each function. 

In this Letter we propose a complementary approach 
to this issue which will finally use only three Re,-- 
independent universal parameters which will be used for 
all the functions discussed. First we ask what could be 
missed in the textbook derivations of the classical log- 
law (1) which may lead to different velocity profile [in- 
cluding possibly the power law (2)]? Our answer is: the 
mean turbulent velocity profile in the entire channel or 
pipe can be described within the traditional approach if 
one realizes how the characteristic length-scale, which has 
physical meaning of the size of energy containing eddies 
£, depends on the position in the flow. Simple scaling 
near the wall, = kz~^, leads to the log-law (1). The 
alternative suggestion of [5], cx (2;+)"(I^6^), leads to 
alternative power-law (2). We see no physical reason why 
£ should behave in either manner. Instead, we propose 
that the eddy size £ should be about z for z <C L, and 
saturate at some level is ^ L approaching the center- 
line, where the effect of the opposite wall is felt. Our 
analysis of DNS data provides a strong support to this 
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FIG. 1: Color online. Left and Right upper panels: comparison of the theoretical mean velocity profiles (red solid lines) 
at different values of Rer with the DNS data for the channel flow [2, 3] (Left panel, grey squares; model with ^buf ~ 49, k = 
0.415, Is = 0.311) and with the experimental Super-Pipe data [4] (middle panel, grey circles; model with ^buf = 46, k = 
0.405, is = 0.275). In orange dashed line we plot the viscous solution = . In green dashed dotted line we present the 
von-Karman log-law. Note that the theoretical predictions with three Rer-independent parameters fits the data throughout 
the channel and pipe, from the viscous scale, through the buffer layer, the log-layer and the wake. For clarity the plots are 
shifted vertically by five units. Lower panel: The Reynolds-stress profiles (solid lines) at Rbt from 394 to 2003 (in channel) 
and from 5050 to 165,000 (in pipe) in comparison with available DNS data (dots) for the channel. 



idea, allowing us to get, within the traditional (second- 
order) closure procedure, a quantitative description of 
the mean shear, S{z) = dV{z)/dz, the kinetic energy 
density (per unit mass), K{z) = /2, and the tan- 

gential Reynolds stress, W{z) = — {uxUz), in the entire 
flow and in a wide region of Rei-, using only three Re,-- 
independent parameters, k, B and (^s ~ 0.311 L for 
the channel and is « 0.275 L for the pipe). 

The closure model should relate three objects: 5+, 
and W'^. The first (exact) relation between these 
objects follows from the Navier-Stokes equation for the 
mean velocity, being the mechanical balance between the 
momentum generated at distance z from the wall, i.e. 
p'{L — z), and the momentum transferred to the wall by 



kinematic viscosity and turbulent transport. In physical 
and wall units it has the form: 

lyS+W = p'{L-z)^S+ + W+ ^ l-C, C = z/L. (3) 

Already in 1877 Boussinesq attempted to close this equa- 
tion by introducing the notion of turbulent viscosity , 
writing W = ^,^3 [10]. Estimating one 
finishes with the closure W'^ = n^l'^^/ K+S'^ . Here i^ 
is a C-dependent characteristic scale of energy contain- 
ing eddies, determining the nonlinear dissipation of W ^ 
and is a constant introduced here for convenience. 
A more careful analysis of the balance equation for W 
(see Ref. [11] and Appendix) that includes the viscous 
dissipation of W , leads to a somewhat more invloved clo- 
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sure for in a form involving an additional universal, 
Rcx-independent dimensionless function of 

r„W+ « K^e+VWS+, r^{z+) ^ (4) 

Here w 49 is a Rcr-indcpendent length that plays a 
role of the crossover scale (in wall units) between the 
buffer and log- law region. In this form, i^{C) oc z 
near the wall, and the choice ~ 0.20 ensures that 

limc^oC(0=C- 

A third relation to supplement Eqs. (3) and (4) is ob- 
tained by balancing the turbulent energy generated by 
the mean flow at a rate SW , and the dissipation at a 
rate £^ = v{\Vu\'^): 

Here the dissipation is estimated via the energy cascade 
over scales involving a characteristic scale of energy con- 
taining eddies, (.^{z) determining the energy transfer 
rate. The constant k,^ will be used to ensure that the 
slope of this function at z+ = is unity. 

Note that in Eqs. (4) and (5) we used a local-balance 
approximation, neglecting the spatial energy flux. This 
approximation is very good in the log-law region but it 
deteriorates near the wall and near the center-line. Nev- 
ertheless for our purposes this has no consequences. Near 
the wall <C 5+ and the local-balance approximation 
plays no role in the exact mechanical balance (3) that 
determines S. For the same reason we also do not need 
to introduce a correction rj^(z+) in Eq. (5) due to the 
direct viscous dissipation (similar to r^{z'^) in Eq. (4) 
since the length scale replacing ^^^j here will be the dis- 
sipative scale ^diss ~ 5 which is entirely buried in the 
region where W and K are small. Near the centerline 
S'^ tends to zero and Eq. (3) determines W'^ « 1 — C, 
which allows an accurate determination of S'^ , because 
we know that and must saturate. 
Profiles of the characteristic length-scales 1^^, 
Now we show that the source of confusion is the assump- 
tion that the relevant length scales can be determined 
a-priori as £~^_ oc {z^)" with a = 1 or a ^ 1. The ac- 
tual dependence £„ and £^ on z and L can be found from 
the data provided by the numerical simulations. Consider 
first £^^-, defined by Eq. (4). We expect that plotting the 
scaling function i?+ /Re,- computed for different values of 
Re^ should collapse the date onto one scaling function. 
The quality of the data collapse for this scaling function 
is presented in Fig. 2, demonstrating the expected satu- 
ration at the center-line. 

The second length-scale, £'^, is determined by the sec- 
ond of Eq. (5). We again expect that i?+/Rer should 
collapse the data obtained from different value of Re,- 
onto one scaling function. In Fig. 2 we demonstrate that 
this scaling function leads to acceptable data collapse 



throughout the channel and for all the four values of Rer 
for which the simulation data are available. 
Solution, Velocity Profiles and Final Scaling 
Function: Solving Eqs. (3) together with S~^W~^ = 

K^^'"^ /{i^K^-t) that follows from Eq. (5), we find 

W+ = {nS+l+fr-J'\ (6) 

where we have defined the von-Karman constant k = 
(k^ Hj^y/'^ « 0.415 and the crucial scaling function ^"'"(C) 
as follows 

^+ = K'iO^tiO]'^' = <lw+'rl/S+ht. (7) 

Note that if one replaces the energy dissipation rate e+ 
by the rate of energy production S'^ and takes r,^, as 
unity this scaling function becomes the Prandtl mixing 
length [1]. However the latter suffers from a non-physical 
divergence at the center-line whereas our length saturates 
to a constant there as it should. 

The convincing data collapse for the resulting function 
£"'"(C)/ReT is shown in Fig. 2, rightmost panel. Substi- 
tuting Eq. (6) in Eq. (3) we find a quadratic equation for 
S with a solution: 

^ _ ^l + {l-0[2n£+{0]yr^{z+fJ^ -I 

^ ~ 2[n£+{C)?/r„{^+Y'^ ' 

To integrate this equation and find the mean velocity 
profile for any value of Re,- we need to determine the 

scaling function £^{C) from the data. A careful analysis 
of the DNS data allows us to find a good one-parameter 
fit for ^+(0 

[-1(^4)]} 

where C = C(l - C/2) and 4 ~ 0.311. The quality of the 
fit is obvious from the continuous line in the rightmost 
panel of Fig. 2. Note that the fit fimction is exactly 
constant at mid channel, with zero slope. This is required 
by symmetry, and will be the reason for our excellent fit 
of data in the wake region. 

Finally the theory for the mean velocity contains three 
parameters, namely 4 together with f^^g (which deter- 
mines B in Eq. (1)) and k. We demonstrate now that 
with these three parameters we can determine the mean 
velocity profile for any value Re^, throughout the chan- 
nel, including the viscous layer, the buffer sub-layer, the 
log-law region and the wake. Examples of the integration 
of Eq. (8) are shown in Fig. 1. It is worthwhile to re- 
iterate that the excellent fits in the viscous and the wake 
regions (superior to the fits presented in [11, 12]), which 
are usually most diflacult to achieve, are obtained here 
due to the correct asymptotics of ^"""(C) at C ^ and 
C, ^ \. In addition, our theory results also in the kinetic 
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FIG. 2: Color online. The scaling function i+{C)/Re-r (Left upper panel), i+iQ/Re^ (Right upper panel) and the final scaling 
function £+{() (L ower panel), as a function of = z/L, for four different values of Re^, computed from the DNS data [2, 3]. 
Note the data collapse everywhere except at 1 where ~ S"*" ^ 1 and accuracy is lost. The green dash line represents 
= (1 — (^/2) with a saturation level 0.5; in orange solid line we show the fitted function Eq. (9) with ^sat ~ 0.311. 



energy, and Reynolds stress profiles which are in a quan- 
titative agreement with the DNS data; for W profiles see 
Fig. 1. 

Conclusions and application to experiments: We 

discussed turbulent channel flow, demonstrating the ex- 
istence and usefulness of a scaling function which 
allows us to get the profiles of the mean velocities for 
all values of Re,- and throughout the channel, in a good 
agreement with DNS. We argued that the controversy be- 
tween power-laws and log-laws is moot, stemming from a 
rough estimate of the scaling function ^^(C)- While this 
function begins near the wall as z+, it saturates later, 
and its full functional dependence on C is crucial for find- 
ing the correct mean velocity profiles. The approach also 



allows us to delineate the accuracy of the log-law presen- 
tation, which depends on and the value of Re^-. For 
asymptotically large Re,- the region of the log-law can 
be very large, but nevertheless it breaks down near the 
mid channel and near the buffer layer, where correction 
to the log-law were presented. 

To show that the present approach is quite general, we 
apply it now to the experimental data that were at the 
center of the controversy [5], i.e. the Princeton Univer- 
sity Superpipe data [4] . In Fig. 1 right panel we show the 
mean velocity profiles as measured in the Superpipe com- 
pared with our prediction using the same scaling function 
£+((■). Note that the data spans values of Re,- from 5050 
to 165000, and the fits with only three Re^-independent 
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constants arc very satisfactory. Note the 2% difFcrence 
in the value of k between the DNS and the experimental 
data; we do not know at this point whether this stems 
from inaccuracies in the DNS or the experimental data, 
or whether turbulent flows in different geometries have 
different values of k. While the latter is theoretically 
questionable, wc cannot exclude this possibility until a 
better understanding of how to compute k from first prin- 
ciples is achieved. 
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the second term is dominant only near the wall where 
^+ = z'^, then r„ ^ 1 + tt^J z'^ . In [12] it was realized 
that this from, which is an interpolation between the near 
wall and the bulk physics, can be modeled in a way that 
reflects better the actual width of the buffer layer, using 
another interpolation formula that reads 



l/n 



(11) 



with n = 2. Best fit to simulational data which is cur- 
rently available is obtained with 5 < n < 7. In this Letter 
we chose n = 6 leading to the second of Eqs. (4). This 
choice simplifies the appearance of the Eqs. (6)- (8). 



Appendix: The exact balance equation for the 
Reynolds shear stress can be found in [1]: P+ + TZ^ 

i Here P+ = -T+yS+ is the production of W+, 
is the redistribution of between other Reynolds 
stress components, £^ is the viscous dissipation of W'^ 
and T+ is the turbulent transport of . Explicit ex- 
pressions for these terms are in [1]. Since Tyy is 0{K), 
we approximate P+ a -K+S+. 11+ = R^'+ + 
[1, 11]. The first term describes the return to isotropy 
, while the second one is responsible for the isotropiza- 
tion of production. A slightly modified Rotta's model 
[13] proposes that i?^,' cx ^/KW/i„. i?'^ is modeled 
according to [1, 14], such that cx K+S+. 

The viscous dissipation e„ = f (dkUxdkUz) is 
0{—i'Wz^'^). As explained in the text, we can neglect 
the non-local term T,^ in the balance for the Reynolds 
stress with impunity. To compensate for its loss in the 
viscous range we increase the estimate {—vWz^"^) by a 
factor Kj K^, , where if* is a dimensional constant [11] 
Eventually, oc —W+VK+/z+'^. Hence, the approxi- 
mate algebraic balance equation for the Reynolds shear 
stress reads: 



-aK+S+ + h- 



■ + CK+S+ « -d- 



(10) 



where a, b,c,d - are positive constants of 0(1). The last 
equation may be rearranged to the form of the fist of Eq. 
(4) but with = 1 + f + , e+/z+^ , e+, = d/b. Since 
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